Search on a Hypercubic Lattice using a Quantum Random Walk: I. d > 2 
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Random walks describe diffusion processes, where movement at every time step is restricted to only 
the neighbouring locations. We construct a quantum random walk algorithm, based on discretisation 
of the Dirac evolution operator inspired by staggered lattice fermions. We use it to investigate the 
spatial search problem, i.e. find a marked vertex on a d-dimensional hypercubic lattice. The 
restriction on movement hardly matters for d > 2, and scaling behaviour close to Grover's optimal 
algorithm (which has no restriction on movement) can be achieved. Using numerical simulations, we 
optimise the proportionality constants of the scaling behaviour, and demonstrate the approach to 
that for Grover's algorithm (equivalent to the mean field theory or the d — > oo limit). In particular, 
the scaling behaviour for d = 3 is only about 25% higher than the optimal d — > oo value. 

PACS numbers: 03.67.Ac, 03.67.Lx 



SPATIAL SEARCH USING THE 
DIRAC OPERATOR 



The spatial search problem is to find a marked object 
from an unsorted database of size TV spread over dis- 
tinct locations. Its characteristic is that one can proceed 
from any location to only its neighbours, while inspect- 
ing the objects. The problem is conventionally repre- 
sented by a graph, with the vertices denoting the loca- 
tions of objects and the edges labeling the connectivity 
of neighbours. Classical algorithms for this problem are 
O(N), since they can do no better than inspect one lo- 
cation after another until reaching the marked object. 
On the other hand, quantum algorithms can do better 
in search problems by working with a superposition of 
states, Grover's algorithm being the prime example 
The spatial search problem has been a focus of inves- 
tigation in recent years using different quantum algo- 
rithmic techniques, both analytical and numerical, and 
in a variety of spatial geometries ranging from a single 
hypercube to regular lattices 0-0] • These studies have 
mainly looked at the asymptotic scaling forms of the al- 
gorithms, and have not varied the database size iV and 
its dimension d independently. In this work, we study 
the specific case of searching for a marked vertex on a 
rf-dimensional hypercubic lattice with N = L d vertices. 
We let N and d be independent, as well as determine the 
scaling prefactors, and thereby develop a broad picture of 
how the dimension of the database (or the connectivity 
of the graph) influences the spatial search problem. 

The quantum algorithmic strategy for spatial search is 
to construct a Hamiltonian evolution, where the kinetic 
part of the Hamiltonian diffuses the amplitude distribu- 
tion all over the lattice and the potential part of the 
Hamiltonian attracts the amplitude distribution toward 
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the marked vertex [8| . The optimization criterion for the 
algorithm is to concentrate the amplitude distribution 
toward the marked vertex as quickly as possible. Grover 
constructed the optimal global diffusion operator, but it 
requires movement from any vertex to any other vertex 
in just one step. That corresponds to a fully connected 
graph, or equivalently the mean field theory limit. When 
movements are restricted to be local (i.e. one can go from 
a vertex to only its neighbours in one step), the search 
algorithm slows down. It is then worthwhile to redesign 
the search algorithm by understanding the extent of slow 
down due to local movements, and how it depends on the 
spatial arrangement of the database. That is the question 
we address quantitatively in this article. As discussed in 
what follows, for the best spatial search algorithms, the 
slow down is in the scaling prefactor for d > 2 and in the 
scaling form for d < 2. 

A diffusion process can be modeled either using con- 
tinuous time and a first order time derivative, or using 
discrete time and a single time step evolution. In its 
discrete version, it is generically described as a random 
walk, whereby a probability distribution evolves in a non- 
deterministic manner at every time step. Such random 
walks have been used to tackle a wide range of graph 
theory problems, usually with a local and translation- 
ally symmetric evolution rule. We use a quantum ver- 
sion of this process, i.e. quantum random walks It 
provides a unitary evolution of the quantum amplitude 
distribution, such that the amplitude at each vertex gets 
redistributed over itself and its neighbours at every time 
step. It is worth noting that quantum random walks are 
deterministic, unlike classical random walks, with quan- 
tum superposition allowing simultaneous exploration of 
multiple possibilities. Several quantum algorithms have 
used them as important ingredients, and an introductory 
overview can be found in Ref . fioj . 

On a periodic lattice with translational symmetry, spa- 
tial propagation modes are characterized by their wave 
vectors k. Quantum diffusion depends on the energy of 
these modes according to U(k,t) = exp(—iE(k)t). The 
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lowest energy mode, k = 0, corresponding to a uni- 
form distribution, is an eigenstate of the diffusion op- 
erator and does not propagate. The slowest propagating 
modes are the ones with the smallest nonzero \k\. The 
commonly used diffusion operator is the Laplacian, with 
E{k) oc \k\ 2 . The alternative Dirac operator, available 
in a quantum setting, provides a faster diffusion of the 
slowest modes with E{k) oc \k\. Its quadratically faster 
spread makes it better suited to quantum spatial search 
algorithms 0, 0] . 

An automatic consequence of the Dirac operator is the 
appearance of an internal degree of freedom correspond- 
ing to spin, whereby the quantum state becomes a multi- 
component spinor. These spinor components can guide 
the diffusion process, and be interpreted as the states of 
a coin [3|, While this is the only possibility for the 
continuum theory, another option is available for a lat- 
tice theory, i.e. the staggered fermion formalism (Tlj . In 
this approach, the spinor degrees of freedom are spread in 
coordinate space over an elementary hypercube, instead 
of being in an internal space. We follow it to construct 
a quantum search algorithm on a hypercubic lattice, re- 
ducing the total Hilbert space dimension by 2 d and elim- 
inating the coin toss instruction. 



A. Bipartite Decomposition 

The free particle Dirac Hamiltonian in d dimensions is 
^froc — —iot ■ V + (3m . (1) 

On a hypercubic lattice, with the simplest discretization 
of the gradient operator, 

V„/(2) = \[f{x + n)-f(x-n)] , (2) 

and a convenient choice of basis, the anticommuting ma- 
trices a, P are spin-diagonalized to location dependent 
signs: 



ifj -> Tip, H -> THT\ T = of/a*- 1 



d ^d-l 
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a n ^]J(-ir, (3^1[(-1T>(3 . (4) 

3=1 3=1 

With this choice, translational invariance holds in steps 
of 2 (instead of 1), and the squared Hamiltonian, 



He, 



-V ■ V + m 2 , 



(5) 



has eigenvalues sin 2 (fc„) + m 2 on a hypercubic lattice. 
Inclusion of a mass term slows down the diffusion process, 
but it also regulates the k = mode. In the present 
article, we set m = 0; in an accompanying article [l2| . 
we investigate the advantage a nonzero mass offers when 
the problem has an infrared divergence. 

Even when the Hamiltonian H is local, the discrete 
time evolution operator U = exp(— iHr) may connect 



arbitrarily separated vertices. To make U also local, we 
split H in to block-diagonal Hermitian parts and then 
exponentiate each part separately For a bipartite 

lattice, partitioning of H in to two parts (which we label 
"odd" and "even" ) is sufficient for this purpose [13] : 

//free = H Q + H e . (6) 

This partition is illustrated in Fig.l for d = 1 and d = 2. 
Each part contains all the vertices, but only half of the 
links attached to each vertex. Each link is associated 
with a term in Hf ree providing propagation along it, and 
appears in only one of the two parts. H thus gets divided 
in to a set of non-overlapping blocks of size 2 d x 2 d (each 
block corresponds to an elementary hypercube on the lat- 
tice) that can be exactly exponentiated. The amplitude 
distribution then evolves according to 



ip(x; t) = W f V>0?; 0) , W = U e U a = e 



-iH e r —iH a T 



(7) 



Each block of the unitary matrices t/ ( e ) mixes ampli- 
tudes of vertices belonging to a single elementary hy- 
percube, and the amplitude distribution spreads in time 
because the two alternating matrices do not commute. 
Note that W does not perform evolution according to 
i/froo exactly. Instead, W = exp(— iH{ lee T) + 0(t 2 ). 
Still, the truncation is such that W is exactly unitary, 
i.e. W = exp(-iHr) for some H. 

We adopt the convention where the walk operator W 
is real. For d = 1, we set 



H a \x) 
HJx) 



(-ir\ x +(-ir) 

(-l)*\ x -(-i)°) 



(8) 
(9) 



The resultant Hamiltonian blocks arc Hf = — 02/2, and 
Hf = CT2/2 when operating on the block with coordi- 
nates flipped in sign. With (H^) 2 = (iff) 2 = \l, 



U ( e ) = cl — 2isH, 



o(e) 



= 1 , 



(10) 



(a) 



(b) 






FIG. 1: Bipartite decomposition of the hypercubic lattice in 
to odd and even parts: (a) for d = 1, (b) for d = 2. 
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where s = sin(r/2) is a parameter that can be tuned. 

These expressions can be rearranged to separate the 
translationally invariant part of the walk from the loca- 
tion dependent signs. Let the amplitude distribution in 
a two spinor (or coin) component notation be 



- ( i>(2x,t) 



\1>(2x-l,t) 



Then the walk operator becomes 

W = UC , C = U* 
U\X) = 



c s 
~s c 



c\X)-sJ2(±a ± )\X±l) 
± 



(11) 

(12) 
(13) 



The operator C = just mixes the components of ^ 
within the block. The first term of the operator U is sta- 
tionary, while the other moves the amplitude to neigh- 
bouring blocks. The movement is accompanied by chi- 
rality flip, a characteristic feature of the Dirac opera- 
tor, because of the raising and lowering Pauli operators 
<t± = (a i ± ia 2 )/2. 

In d dimensions, the preceding structure generalises 
to 2 d x 2 d Hamiltonian blocks that can be expressed in 
terms of tensor products of Pauli matrices [151 ] . When 
operating on hypercubes with coordinate labels {0, l}® d , 



H, 



1 d 



<x 2 



<8>0-i) 



(14) 



and Hg = —H^ when operating on a hypercube with 
all coordinates flipped in sign. Note that Eq. (fT4"|) is the 
sum of d terms, which mutually anticommute and whose 
squares are proportional to identity. This is the charac- 
teristic signature of the Dirac Hamiltonian, Eq.([T|). The 
block-diagonal matrices satisfy H% 
ponentiate to 



H 



\I, and ex- 



U, 



o(e) 



cl — is 



Vd 



H. 



o(e) 



(15) 



The parameter s — sin(v / dr/2) can be optimised, to 
achieve the fastest diffusion across the lattice. Note that 
the d-dimensional walk is not the tensor production of d 
one-dimensional walks. 

Again separating the degrees of freedom within each 
local hypercube as x = X®{0,-l}® d , the location depen- 
dent signs can be rewritten as operators acting on the 
spinor (or coin) degrees of freedom. The walk operator 
takes the form: W = UC, C = U* , and 



U\X) = c\X) 



(16) 



• 5353(i«(*-i) ® (± CT± ) *f d -»)\x ± j) 

± j=l 



B. Search Algorithm 

To search for a marked vertex, say the origin, we at- 
tract the quantum random walk toward it by adding a 
potential to the free Hamiltonian, 



y _ f j3Vo 5g >0 (gravitational) , 
\ Vo Ss,o (electromagnetic) . 



(17) 



Exponentiation of this potential produces a phase change 
for the amplitude at the marked vertex. For the steep- 
est attraction of the quantum random walk toward the 
marked vertex, it is optimal to choose Vq so as to make 
the corresponding evolution phase maximally different 
from 1, i.e. e~ lV ° T = — 1. The sign provided by (3 does 
not matter in this case, and the evolution phase becomes 
the binary oracle, 



R = I-2 



(18) 



The search algorithm alternates between diffusion and 
oracle operators — a discrete version of Trotter's formula 
involving H , H e and V — yielding the evolution 



i)(x;t u t 2 ) = [W* 1 R] t2 ip(x; 0, 0) . 



(19) 



Here ti is the number of search oracle queries, and t\ is 
the number of walk steps between queries. Both should 
be minimised to find the the quickest solution to the 
search problem. 

Since this algorithm iterates a unitary operator on the 
initial state, it produces periodic results. As a function of 
the number of iterations, the probability of being at the 
marked vertex first increases, reaches a peak value P after 
ti search oracle queries, then decreases toward zero, and 
thereafter keeps on repeating this cycle. Grover's algo- 
rithm is designed to evolve the quantum state in a two- 
dimensional Hilbert space (spanned by the initial and 
marked states), and is able to reach P = 1. That does 
not hold for spatial search, and the resultant values of P 
are less than 1. The remedy is to augment the algorithm 
by the amplitude amplification procedure [l6j . to find 
the marked vertex with 0(1) probability. The complex- 
ity of the algorithm is then characterised by the effective 
number of search oracle queries, t^/yP- 



C. Complexity Bounds 

Spatial search obeys two simple lower bounds. One 
arises from the fact that while the marked vertex could 
be anywhere on the lattice, one step of the local walk 
can move from any vertex to only its nearest neighbours. 
Since the marked vertex cannot be located without reach- 
ing it, the worst case scenario on a hypercubic lattice re- 
quires Q,(dL) steps. This bound is the strongest in one 
dimension, where quantum search is unable to improve 
on the O(N) classical search. 

The other bound follows from the fact that spatial 
search cannot outperform Grover's optimal algorithm, 
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and so it must require il(y/N) oracle queries. That is 
independent of d, and stronger than the first bound for 
d > 2. As a matter of fact, the first bound weakens with 
increasing d, as the number of neighbours of a vertex in- 
crease and steps required to go from one corner of the 
lattice to another decrease. In analogy with mean field 
theory behaviour of problems in statistical mechanics, 
we can thus expect that in the d — ¥ oo limit the local- 
ity restriction on the walk would become irrelevant and 
the complexity of spatial search would approach that of 
Grover's algorithm. (Note that the maximum value of d 
is log 2 N for finite N.) Our results in this work demon- 
strate that not only is the Q(y/N) complexity scaling 
achievable for d > 2, but one can get pretty close to the 
asymptotic prefactor ir/4 as well. 

Combined together, the two bounds make the complex- 
ity of spatial search fl(dN 1 ^ d : y/~N). The two bounds are 
of the same magnitude, Q(yN), for the critical case of 
d = 2. It is a familiar occurrence in statistical mechan- 
ics of critical phenomena that an interplay of different 
dynamical features produces logarithmic correction fac- 
tors in critical dimensions due to infrared divergences. 
The same seems to be true for spatial search in d = 2, 
where algorithms are slowed down by extra logarithmic 
factors [3-l5(. It is an open question to design algorithms, 
perhaps using additional parameters, that suppress the 
logarithmic factors as much as possible. We look at the 
special d = 2 case in an accompanying article. 

The way the lower bounds arise in spatial search also 
illustrates an interplay of two distinct physical princi- 
ples, special relativity (or no faster than light signaling) 
and unitarity. It is well-known that the two are compati- 
ble, but just barely so, in the sense that physical theories 
that are unitary but not relativistic or relativistic but not 
unitary exist. Furthermore, the best algorithms should 
arise in a framework that respects both the principles, 
i.e. relativistic quantum mechanics. In spatial search, 
locality of the quantum walk and the square-root speed 
up produced by the Dirac equation (through change in 
the dispersion relation) follow from special relativity. On 
the other hand, optimality of square-root speed up pro- 
vided by Grover's algorithm is a consequence of unitarity 
[l7l [l8j . Why these two distinct physical principles lead 
to the same scaling constraint is not understood, and is 
a really interesting question to explore. At present, how- 
ever, we have nothing more to add. 

We point out that preparation of the unbiased initial 
state for the search problem, i.e. the translationally in- 
variant uniform superposition state \s) = ^2 X \x)/y/N, is 
not at all difficult using local directed walk steps. For 
instance, one can start at the origin, step by step trans- 
fer the amplitude to the next vertex along an axis, and 
achieve an amplitude L~ x / 2 at all the vertices on the axis 
after L steps. Repeating the procedure for each remain- 
ing coordinate direction produces the state \s) after dL 
steps in total. Clearly, this preparation does not add to 
the complexity of the algorithm. 



II. OPTIMISATION OF PARAMETERS 

The fastest search amounts to finding the shortest uni- 
tary evolution path between the initial state \s) and the 
marked state |0). This path is a geodesic arc on the uni- 
tary sphere from \s) to |0), and Grover's optimal algo- 
rithm strides along this path perfectly. In our algorithm, 
Ea.(fT9)). we have replaced Grover's optimal diffusion op- 
erator, G = 2\s)(s\ — 1, with the walk operator W* 1 . So 
to optimise our algorithm, we need to tune the param- 
eters appearing in W tx , i.e. s (or c) and ti, such that 
W ll R approximates a rotation in the two-dimensional 
|s)-|0) subspace by the largest possible angle. 

The operator G has |s) as an eigenstate with eigen- 
value 1. All other states orthogonal to \s), in particular 
the combination \s±) = {\s)-y/N\6))/y/N - 1 in the \s)- 
|0) subspace, are its eigenstates with eigenvalue —1. The 
walk operator W = U e U also has \s) as an eigenstate 
with eigenvalue 1. However, its other eigenvalues are 
spread around the unit circle, and evolution under W tx R 
does not remain fully confined to the |s)-|0) subspace. If 
the outward diffusion is not controlled, the probability of 
remaining in the |s)-|0) subspace would decay exponen- 
tially with the number of iterations. To reach the marked 
state with a constant probability, therefore, parameters 
must be tuned so that part of the outward diffusion sub- 
sequently returns to the |s)-|0) subspace. 

A. Analytical Criteria 

The optimisation condition that makes W* 1 as close 
an approximation to G as possible can be formulated in 
several different ways. Though they all lead to the same 
result, describing them separately is instructive. 

(1) Maximise the overlap of W' 1 and G, i.e. 

Tr{W ll G) = Tr(2\s)(s\-W tl ) = 2-^{i\W tl \i) . (20) 

i 

(i\W l1 \i) is the amplitude for the random walk to return 
to the starting state \i) after t\ time steps. Since W is 
translationally invariant along each coordinate direction 
in steps of 2, we need to evaluate (ilVF* 1 \i) only for states 
|z) in an elementary hypercube. Now, any closed path 
on a hypercubic lattice takes an even number of steps 
along each coordinate direction, and all location depen- 
dent propagation signs (cf. Eq.(|l])) get squared to 1 in 
the process. As a result, (zlFK* 1 |i) is independent of \i), 
and the optimisation condition becomes the minimisation 
of A(ti) = (6\W tl \d). 

(2) Make \s±) approximate an eigenvector of W fl , with 
an eigenvalue approaching —1, i.e. minimise 

(sx\W t i\s x ) = w ?- Y (-l + NA(t 1 )) , (21) 

which happens to be the same condition as given earlier. 
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In addition, unitarity of the evolution operator implies exp(±ir Jj^j sin 2 kj), with the eigenphases in the inter 



\{ S ±\W^\ S1 _)\<1 



-l + jj<A(h)<l. (22) 



(3) Maximise the rotation provided by the operator 
W^R along the |s)-|0) subspace. This rotation is given 
by the projection 



Proj 




The factor on the right in Eq.([23j) is precisely the ro- 
tation matrix produced by the operator GR in Grover's 
algorithm, i.e. rotation by an angle 2 sin -1 (1 / '\N) in the 
|s)-|0) subspace. So the total projection can be viewed 
as first applying an iteration of Grover's algorithm to the 
state, and then shrinking the |sj_) component by a fac- 
tor (1 — N A(ti)) I (N — 1). The shrinking disappears and 
the projection tends to an exact orthogonal transforma- 
tion, as A(ti) approaches its lower bound (2 — N)/N. It 
also follows that when the shrinking is controlled, cutting 
down outward diffusion from the |s)-|0) subspace, the al- 
gorithm will have the same 0(y/~N) scaling as in case of 
Grover's algorithm. 

The minimisation condition for A{t{) is obtained here 
only using the translational invariance of the walk op- 
erator W, together with W\s) = \s) and (s\W = (s\. 
It is easy to perform this optimisation numerically, as 
a function of the parameters s and ti, and our results 
are described later. Before that, an understanding of the 
eigenspectrum of W gives us some idea regarding what 
to expect. 

Since iJ 2 = -ff 2 = 4l, eigenvalues of H a and H e 
are ±y/d/2. Moreover, H and H e are traceless, so 
that half the eigenvalues are yd/2 and the other half 
— yd/2. Consequently, eigenvalues of U and U e are 
c± is = exp(±z-\/rfT /2), equally divided among complex 
conjugate pairs. Furthermore, H r e \ and U / e \ have the 
same eigenvectors. 

Since U and U e do not commute, we do not have a 
general solution for the eigenspectrum of W. Still it is 
straightforward to work it out in the continuum time 
limit, t — Y 0. i/free is diagonalised by a Fourier transform, 



i d 



\x + 2h) - 2\x) + \x-2h) 



]Tsin 2 (fc n )|fc) , 



(24) 



71=1 



with a uniform distribution in the Brillouin zone. Up to 
0(t 2 ), therefore, the eigenvalues of W(~ e~ lHt " cT ) are 



val [-rVd^rVd]. The average value of Hf ree is d/2, and 
the corresponding eigenvalues of W are exp(±zr-y/ d/2). 
We find that this average value plays an important role 
in determination of the optimal parameters. 



B. Numerical Optimisation 

The considerations of the previous subsection focused 
on what happens during a single iteration of the algo- 
rithm. Over many iterations, there is outward diffusion 
of the amplitude distribution from the |s)-|0) subspace, 
and subsequent partial return. That plays an important 
role in the overall success of the algorithm, but we have 
not been able to estimate it analytically. To study the 
performance of the whole algorithm, therefore, we nu- 
merically tune the parameters s and t\ to (a) maximise 
the probability of reaching the marked vertex, and (b) 
minimise the number of search oracle queries required to 
reach that stage. The ensuing results of our optimisation 
of P, t-x and A(t\) are described in what follows. 

We first looked at how the probability distribution 
evolves as a function of time during the search process. 
To illustrate our results, we take s = c = l/v2 that 
corresponds to maximal mixing between even and odd 
vertices, and analyse the evolution as a function of t\ 
and t% on a 32 3 lattice. As displayed in Fig. 2, the proba- 
bility at the marked vertex undergoes a cyclic evolution 
pattern as a function of time t?.- We observe that the 
evolution is essentially sinusoidal for t\ = 1,2,3, and 
persists for more than 30 cycles without any visible devi- 
ation. This behaviour implies that a constant fraction of 
the quantum state remains in the |s)-|0) subspace, and 
the operator W ll R rotates it at a uniform rate within the 
subspace, fully matching the exact solution of Grover's 
algorithm. The difference among t\ — 1,2,3 is that the 
cycle time and the peak probability vary, meaning that 
the fraction of the quantum state present in the |s)-|0) 
subspace depends on t\. 

For ti > 3, we find that the evolution rapidly loses 
its sinusoidal nature and the peak probability plummets. 
The change in the pattern is so drastic that the algorithm 
no longer performs a reasonable search. This behaviour 
implies that most of the quantum state has moved out 
of the |s)-|0) subspace and only a negligible component 
is left behind. An appropriate choice of t± is thus crucial 
for constructing a successful search algorithm. 

We also show in Fig. 3 what the probability distribu- 
tion looks like at the instance when the probability at 
the marked vertex attains its peak value. We note that 
the distribution is strongly peaked around the marked 
vertex in an essentially uniform background. The peak 
is sharper when the peak probability is larger, and it al- 
most disappears for t\ > 3. This pattern suggests that 
the part of the quantum state that diffuses outside the 
|s)-|0) subspace is almost featureless. 
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Combining our observations in this particular example, 
we find that the best results — the shortest period, the 
largest peak probability and the sharpest peak — clearly 
correspond to t\ = 3. The results improve as t\ increases 
from 1 to 3, and rapidly deteriorate thereafter. The im- 
portant fact that all optimal features appear simultane- 
ously at t\ = 3 simplifies tuning of the parameters. 

Next we scanned through values of s, for different val- 
ues of ti and in different dimensions, to determine the 
optimal values at which the peak probability P at the 
marked vertex is the largest and the time £2 required to 
reach it is the smallest. A typical situation, with t\ = 3 
on a 32 3 lattice, is shown in Fig. 4. We notice that con- 
tinuous P has a smoother behaviour than discrete £2, 
and maximisation of P and minimisation of £2 occur at 
roughly the same value of s. So we simplify matters, and 
select maximisation of P as our optimality condition. 




500 1000 1500 

Time 






FIG. 2: (Color online) Time evolution of the probability at the 
marked vertex for the spatial search problem on a 32 3 lattice, 
s = l/\/2 and ti increases from 1 (top) to 4 (bottom). 



FIG. 3: (Color online) Probability distribution for the spatial 
search problem on a 32 3 lattice, at the instance when the 
probability at the marked vertex (16, 16, 16) attains its peak 
value. The displayed results are for the z = 16 slice, s — 1/V2 
and t\ increases from 1 (top) to 4 (bottom). 
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A comparison of our scans, for t\ = 3 and different 
number of dimensions, is presented in Fig. 5. We find 
that at fixed t%, the shape of the P vs. s curve is al- 
most dimension independent, and the optimal value of s 
is approximately the same for all dimensions. We also 
observe that for fixed s, P roughly scales as 2~ d . In the 
staggered fermion implementation of the Dirac operator, 
different vertices of an elementary hypercube correspond 
to different degrees of freedom [11] . So our observation 
suggests that to a large extent only the degree of free- 
dom corresponding to the marked vertex evolves during 
the search, while other degrees of freedom remain spec- 
tators. Moreover, the feature that the maximum value 
of P approaches 2™ d with increasing d implies that the 
efficiency of our algorithm increases with d. 

All our quantitative results for optimal values of s are 
collected in Table I. We note that the best performance 
parameters, P and £2, depend on d but hardly on t\ or 
s. (With increasing ti, there is a slight increase in P and 
a small decrease in t-i , but that is a rather tiny improve- 
ment.) More interestingly, we find a relation between the 
optimal values of s and t\, in terms of the variable 

9 = \/2ti sin" 1 s = hr^/dpl . (25) 

Since sin -1 s describes the unitary rotation performed 
by the single step operators U r e ), and t\ is the number 
of walk steps, 9 is a measure of the unitary rotation per- 
formed by the operator IF* 1 . The corresponding operator 
in Grover's optimal algorithm, G, is a reflection providing 
a phase change of tt. We find that our optimal W fl is also 
close to a reflection, in a sense, with 9 ss tt. Unlike the 
discrete ±1 eigenvalues of G though, eigenvalues of IF* 1 
form a distribution. We do not know this distribution in 
general, but we know it in the limit r — ¥ 0, as described 
at the end of Subsection II. A. This limit corresponds to 
small s and large t\, and Table I shows that 9 indeed 
gets close to tt in this limit. To be precise, this normal- 
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FIG. 4: (Color online) Optimisation of s for t\ = 3 on a 32 3 
lattice, using peak probability P as well as oracle queries t 2 . 



TABLE I: Results of optimal parameter determination for 
different t\ and in different dimensions. 
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0.0284 


3.242 
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5 


4 


0.5376 
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0.0286 


3.211 


0.5292 


-0.8636 


3.155 






8 


0.2726 


147 


0.0288 


3.124 


0.2710 


-0.8718 


3.105 






20 


0.1108 


147 


0.0288 


3.140 


0.1092 


-0.8739 


3.095 


6 


8 


3 


0.6891 


51 


0.0145 


3.225 


0.6913 


-0.8778 


3.238 




20 


0.1106 


51 


0.0147 


3.135 


0.1094 


-0.8951 


3.101 


7 


8 


3 


0.6932 
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0.0073 


3.250 


0.6937 


-0.8949 


3.252 




20 


0.1097 


102 


0.0074 


3.109 


0.1098 


-0.9102 


3.112 



isation of 9 implies that the optimal W fl = e~' lHrtl is a 
reflection operator for the "average eigenvalue" of the as- 
sociated H 2 . Appearance of the average eigenvalue in the 
normalisation indicates that all spatial modes contribute 
to the search process with roughly equal strength. 

Table I also lists our results for the variable A(t\) de- 
fined in Subsection II. A, which characterises how close 
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s 

FIG. 5: (Color online) Optimisation of s for t\ — 3 in different 
dimensions. Lattice sizes L — 32, 16, 16, 8, 8 were used for 
d = 3, 4, 5, 6, 7, respectively. 
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W fl is to G. We computed A(ti), by starting with unit 
amplitude at the marked vertex, and then measuring the 
amplitude at the marked vertex after t\ walk steps. At 
the beginning of this subsection we mentioned that check- 
ing A(t±) is not as elaborate an optimisation test as us- 
ing the full search algorithm. Nevertheless, the results 
corroborate our optimisation of spatial search: The min- 
imum value of A(ti) is reasonably close to —1, and occurs 
at approximately the same value of s (and hence 9) that 
provides the optimal search (the mismatch decreases with 
increasing d). We also observe that A(t{) gets closer to 
— 1, suggesting that the algorithm becomes more efficient, 
with increasing d as well as t%. 

Given the strong correlation between the optimal val- 
ues of s and t\ , the computational cost of the algorithm 
is minimised by choosing the smallest t%. For t\ = 1, 
Eq. ([25|) does not have a solution with 8 ps tt, the optimal 
value of s is very close to 1, and our algorithm does not 
work well there. Then t\ = 2 is the choice that min- 
imises the number of walk steps in the algorithm. Even 
when the emphasis is on optimising P and t2, most of 
the improvement can be obtained by going from t\ = 2 
to t\ = 3, and we need not go beyond that. 



III. SIMULATION RESULTS 

We carried out numerical simulations of the quantum 
spatial search problem, with a single marked vertex and 
lattice dimension ranging from 3 to 9. The matrices Z7 ( e ) 
are real with our conventions, and that was convenient 
for simulations. The optimal values of s-t% pairs depend 
a little on the lattice size and dimension. We are only 
interested in the asymptotic behaviour of the algorithm, 
however, so we used the same parameter values for all 
lattice sizes and dimensions: s — 0.9539 for t\ = 2, s = 



1/V2 for ti = 3 and s = 0.5410 for t x = 4. 

To extract the scaling properties of the algorithm, we 
studied the behaviour of P and t 2 , as a function of L and 
d. We found that the finite size effects in our data, as 
L — > oo at fixed d, are well described by series in inverse 
powers of L. In Figs. 6-7, we show some of our results 
with the fitting functions: 



P = ai 



h 

L 



N 



a 2 



L 



(26) 



All the fit parameters are listed in Table II, where error 
refers to the r.m.s. deviation of the data from the fit. We 
do not have sufficient data for d = 8 and d = 9, to make 
accurate fits or to see the asymptotic behaviour. Our 
analysis, therefore, relies on the patterns seen for d — 3 
to d = 7, while we use d = 8 and d — 9 values only for 
consistency checks. 

We observe the following features: 

(1) There is not much difference among the a\ and a 2 
values for different optimised s-t\ pairs. This confirms 
our earlier inference that implementing our algorithm for 
ti = 2 minimises the running cost. 

(2) P approaches a constant as L — > oo, with a\ decreas- 
ing approximately as 2~ d with increasing d. As men- 
tioned before, this is a consequence of the fact that in 
the staggered fermion implementation of the Dirac op- 
erator, only the degree of freedom corresponding to the 
location of the marked vertex on the elementary hyper- 
cube participates in the search process. Since the number 
of vertices in an elementary hypercube is 2 d , P is upper 
bounded by 2~ d . Our algorithm achieves that bound up 
to a constant multiple close to 1. 

(3) b\ decreases toward zero with increasing d, which 
is also true for I&2I to a large extent. Thus the small- 
est d — 3 has the largest finite size corrections to the 
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FIG. 6: (Color online) Peak probability at the marked vertex 
as a function of database size (ranging from 2 8 to 2 27 ) in differ- 
ent dimensions. The points are the data from the simulations 
with ti = 3, and the curves are the fits P — a\ + (bi/L). 



FIG. 7: (Color online) Number of search oracle queries as a 
function of database size (ranging from 2 s to 2 27 ) in different 
dimensions. The points are the data from the simulations 
with ti = 3, and the curves are the fits fa/y/N) = a2 + (&2/£). 
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TABLE II: Fit parameters for peak probability P and search oracle queries tijvN vs. L. 



s 


ti 


d 


L 


Ol 


h 


Error 


0,2 


b 2 


Error 




0.9539 


2 


3 


64,128,256,512 


0.0911 


0.0964 


2.43X10" 5 


0.3237 


0.1440 


4.57xl0~ 4 


1.072 


4 


16,32,48,64 


0.0522 


0.0100 


3.02X10" 5 


0.2167 


-0.0832 


1.05 xl(T 3 


0.948 


5 


12,16,24,32 


0.0275 


0.0009 


1.70X10" 6 


0.1486 


-0.0243 


2.39xl0 -4 


0.896 


6 


12,14,16,18 


0.0141 


0.0001 


2.31X10" 7 


0.1043 


-0.0208 


1.78X10 -4 


0.878 


7 


6,8,10 


0.0072 


-0.0004 


1.40xl0 -6 


0.0757 


-0.0338 
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6,8 
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1.23xl0 -3 


1.010 


4 


16,32,48,64 


0.0542 


0.0076 


1.46xl0" 5 
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0.0283 
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0.1470 


-0.0300 


2.12xl0~ 4 
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6 


12,14,16,18 


0.0145 


0.0001 


5.79xl0" 7 


0.1035 


-0.0269 


1.88xl0~ 4 


0.860 


7 


6,8,10 


0.0074 


-0.0003 


1.46xl0" 6 


0.0750 


-0.0296 


3.39xl0~ 4 


0.872 


8 


6,8 


0.0037 




3.00xl0" 6 


0.0498 




4.55xl0~ 5 


0.819 


9 
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0.0019 






0.0353 






0.810 


0.5410 


4 
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64,128,256,512 
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0.0936 


1.84xl0" 5 


0.3123 
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8.46 xl0~ 4 
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2.24 xl0~ 5 


0.2103 


-0.0500 


3.57xl0 -4 


0.898 


5 


12,16,24,32 


0.0285 


0.0013 


8.05xl0~ 6 


0.1455 


-0.0142 


2.25xl0~ 4 


0.862 


6 


12,14,16,18 


0.0146 


0.0002 


7.67xl0" 7 


0.1015 


0.0043 


6.28 xlO -5 


0.840 


7 


6,8,10 


0.0074 


-0.0003 


2.86xl0" 6 


0.0733 


-0.0194 


2.01xl0~ 4 


0.852 


8 


6,8 


0.0037 




4.50xl0" 6 


0.0505 




4.39xl0~ 4 


0.830 


9 


6 


0.0019 






0.0353 
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asymptotic L — > oo behaviour. This suggests that as the 
number of neighbours of a vertex increase, the finite size 
of the lattice becomes less and less relevant. 
(4) t 2 is proportional to ViV as L — > oo, with a 2 de- 
creasing approximately as 2~ d l 2 with increasing d. Thus 
t 2 scales with the database size as (L/2) d / 2 = \JN/2 d . 
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FIG. 8: (Color online) Behaviour of the asymptotic fit param- 
eters as a function of dimension. The data points are from 
Table II, with ti — 3, and the curves are the fits described in 

Eas. ipnm 



This is in agreement with the fact that, with only one 
degree of freedom per elementary hypercube participat- 
ing in the search process, the effective number of vertices 
being searched is AT/2 . 

(5) In the combination t^/yP, describing the effective 
number of oracle queries for our algorithm, the factors of 
2 d arising from the number of vertices in an elementary 
hypercube fully cancel. That makes t2/VP scale as y/N, 
and (Z2 / y/ai describes the complexity scaling of our al- 
gorithm. As listed in the rightmost column of Table II, 
o-i 1 1 \fo-i shows a gradual decrease with increasing d, and 
our results are just above the corresponding value 7r/4 for 
Grover's optimal algorithm. More explicitly, our result 
for d = 3 is about 25% above 7r/4, decreasing to about 
10% above it/ 4 for d = 7. These values imply that our 
algorithm improves in efficiency with increasing d, and 
is not too far from the optimal behaviour even for the 
smallest dimension d — 3. 

To analyse the dimension dependence of our algorithm 
in more detail, we fit the values in Table II to the func- 
tional forms, 



log 2 a\ = ci + did , log 2 a 2 = c 2 + d 2 d 



02_ 

/oi 



C3 



d 



(27) 



(28) 



The results for t\ — 3 are displayed in Fig. 8, showing 
that the fits work reasonably well. 
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TABLE III: Fit parameters for log 2 ai, log 2 a 2 and a 2 / 'yfai vs. d. 



s 


ti 


Fitted d 


Cl 


di 


Error 


C2 


d 2 


Error 


C3 


d 3 


Error 


0.9539 


2 




-0.553 


-0.938 


4.77X10" 2 


-0.085 


-0.526 


2.53X10" 2 


0.721 


0.995 


1.71 xlO" 2 


l/v/2 


3 


3,4,5,6,7,8 


-0.458 


-0.947 


4.48X10" 2 


-0.138 


-0.521 


2.57xl0~ 2 


0.729 


0.791 


1.63xl0~ 2 


0.5410 


4 




-0.421 


-0.951 


3.78X10 -2 


-0.151 


-0.521 


2.09X10" 2 


0.725 


0.761 


1.32X10" 2 



The fit parameters for the three different optimised s- 
ti pairs are listed in Table III. They support our earlier 
scaling observations: d\ close to —1 means that P scales 
roughly as 2~ d , and d 2 close to —0.5 means that t 2 scales 
roughly as \J N/2 d . Interpreting Grover's algorithm as 
the d — ¥ 00 limit of the spatial search algorithm, we ex- 
pect C3 to be close to 7r/4. It turns out to be somewhat 
smaller, but is not very accurately determined due to siz- 
able correction from d$ for the values of d that we have. 
We believe that simulations for larger lattice sizes and 
higher dimensions can remedy this problem. 

To look at our data in a different manner, we analyse 
the scaling behaviour of t 2 /\/NP as a function of the 
database size. In Fi g. 9, we plot a subset of the data, il- 
lustrating how t^j v NP behaves for fixed lattice size L 
when the database size is changed by changing d. The ap- 
proach to the asymptotic behaviour is particularly clear 
in this form, with t 2 /\/ NP decreasing toward 7r/4 as the 
database size increases. It is also obvious that, for fixed 
N, our algorithm improves in efficiency with decreasing 
L (or equivalently with increasing d). This is a conse- 
quence of the increase in the number of neighbours of 
a vertex with increasing d. An important implication is 
that, given a database of size N, it is best to implement 
the algorithm using the smallest L (and hence the largest 
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FIG. 9: (Color online) t 2 /\/ NP as a function of the database 
size for different lattice sizes. The points are the data for 
t± = 3, and the curves are the fits fa/vWP) = ai + (bi/d). 
Also shown is the limiting value corresponding to Grover's 
optimal algorithm. 



d) available. For our staggered fermion implementation 
of the Dirac operator, L has to be even, and the small- 
est L possible is 4. Because of small values of d that we 
have worked with, and discrete nature of t 2 , our data for 
L = 4 do not show smooth asymptotic behaviour. We 
therefore omitted L = 4 from our analysis, and our best 
results are for L — 6. 

We can fit the data reasonably well as 



= ai + — 



(29) 



The fit parameters for the three different optimised s-ti 
pairs are listed in Table IV. We note that a; is close to 
7r/4 when data for d > 6 arc available, but it cannot be 
accurately determined using data for small d only due to 
sizable corrections from bi. 



IV. RESULTS FOR MULTIPLE 
MARKED VERTICES 

We also carried out a few simulations of spatial search 
with more than one marked vertex. The change from 
the single marked vertex algorithm is that the search or- 
acle now flips the sign of amplitudes at several distinct 
vertices, say M. To see clear patterns in the results, we 
needed to keep the marked vertices well separated and 
avoid interference effects. 

With the staggered fermion implementation of the 
Dirac operator, it matters whether the marked vertices 



TABLE IV: Fit parameters for t 2 /VNP vs. d at fixed L. 



s 


h 


L 


Fitted d 


ai 


bi 


Error 






6 


5,6,7,8,9 


0.766 


0.551 


1.06x10" 


2 


0.9539 


2 


8 


4,5,6,7,8 


0.805 


0.293 


5.31x10" 


4 






16 


3,4,5,6 


0.663 


1.132 


1.73x10" 


2 






32 


3,4,5 


0.625 


1.304 


5.77x10" 


3 






6 


4,5,6,7,8,9 


0.819 


0.006 


1.28x10" 


2 


1/V2 


3 


8 


4,5,6,7,8 


0.803 


0.234 


3.28x10" 


3 




16 


3,4,5,6 


0.773 


0.471 


6.46x10" 


3 






32 


3,4,5 


0.724 


0.705 


3.23x10" 


3 






6 


4,5,6,7,8,9 


0.837 


-0.092 


1.35x10" 


2 


0.5410 


4 


8 


4,5,6,7,8 


0.791 


0.272 


3.79x10" 


3 






16 


3,4,5,6 


0.767 


0.450 


1.29x10" 


3 






32 


3,4,5 


0.711 


0.726 


1.13x10" 


3 
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TABLE V: Spatial search results for multiple marked vertices in d — 3. 



s 


tl 


M 


Marked vertices 


P 


ti 






1 


(32,32,32) 


0.09829 


161 






2 


(0,32,32), (32,32,32) 


0.04919, 0.04919 


112, 112 


1/V2 


3 




(0,32,33), (32,32,32) 


0.09868, 0.09790 


161, 161 








(0,0,0), (16,16,16), (32,32,32) 


0.03530, 0.03264, 0.03082 


94, 92, 92 






3 


(0,0,1), (16,16,16), (32,32,32) 


0.09380, 0.05590, 0.04507 


161, 117, 109 








(0,0,1), (16,16,16), (33,32,32) 


0.09298, 0.09838, 0.10347 


157, 161, 162 



belong to the same corner of the elementary hypercube 
or to different corners. So we considered various hyper- 
cube locations for the marked vertices. In Table V, we 
compare our illustrative results for two and three marked 
vertices with those for a single marked vertex. The val- 
ues for peak probability P and corresponding time £2 
were obtained with s = 1/V2 and ti — 3 on a 64 3 lattice. 

When the marked vertices belong to the same corner of 
the elementary hypercube (e.g. (0,32,32) and (32,32,32)), 
they involve the same degree of freedom. In that case, 
we observe that P decreases by a factor of M and t% de- 
creases by a factor of On the other hand, when 
the marked vertices belong to different corners of the ele- 
mentary hypercube (e.g. (0,32,33) and (32,32,32)), they 
involve different degrees of freedom. In that case, we 
find that P and £2 are more or less the same as for the 
single marked vertex case. This pattern is also followed 
when some marked vertices are at the same corner of the 
elementary hypercube and some at different corners. 

We can infer the following equipartition rule from these 
observations. A fixed peak probability (upper bounded 
by 2~ d ) is available for each corner degree of freedom 
of the elementary hypercube. Multiple marked vertices 
corresponding to the same degree of freedom share this 
peak probability equally, while different degrees of free- 
dom evolve independently of each other. Apart from this 
caveat, the dependence of £2 on M for our spatial search 
algorithm is the same as in the case of Grover's algorithm, 



i.e. proportional to y/N/M. 



V. SUMMARY 

We have formulated the spatial search problem us- 
ing staggered fermion discretisation of the Dirac oper- 
ator on a hypercubic lattice in arbitrary dimension. Our 
search strategy, Ea. (fT9|) . mimics that of Grover's algo- 
rithm. The locality restriction on the walk operator be- 
comes less and less relevant as the dimensionality of the 
space increases, and we expect to reach the optimal scal- 
ing behaviour of Grover's algorithm as d — > 00. 

We have verified our theoretical expectations by nu- 
merical simulations of the algorithm for d = 3 to d = 9. 
We clearly demonstrate the approach to the d — » 00 
limit, and find that even for d = 3 our algorithm is only 
25% less efficient. From a computational cost point of 
view, for a fixed database size N, our best parameters 
are t% = 2 and the largest d available. With the stag- 
gered fermion implementation, our algorithm works in 
total Hilbert space dimension N . That is a big reduction 
from the total Hilbert dimension 2 d N required by the 
algorithms that use coin or internal degrees of freedom 
(for the best d = 0(lniV), 2 d = O(N) is substantial), 
and makes our exercise worthwhile. 
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